/*
  Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of the static members of the grain_size classes.

// $Id$

#include "grain_size.h"

std::map<std::string, mcrx::wd01_graphite_distribution::params> mcrx::wd01_graphite_distribution::paramsets;
std::map<std::string, mcrx::wd01_silicate_distribution::params> mcrx::wd01_silicate_distribution::paramsets;

void mcrx::wd01_graphite_distribution::set_params(const params& p) 
{
  alpha_ = p[0];
  beta_ = p[1];
  a_t_ = p[2];
  a_c_ = p[3];
  C_ = p[4];
  b_C1_ = p[5] * 0.75;
  b_C2_ = p[5] * 0.25;
  a0_1_ = p[6];
  a0_2_ = p[7];
  sigma_1_ = p[8];
  sigma_2_ = p[9];
};


mcrx::wd01_graphite_distribution::
wd01_graphite_distribution (const std::string& model,
			    const T_unit_map& units) : 
  wd01_distribution(units) 
{
  if (paramsets.empty()) {
    // initialize the paramsets
    // alpha, beta, a_t, a_c, C, b_C, 
    // a0_1, a0_2, sigma_1, sigma_2
    params p( -1.54, -0.165, 0.0107e-6, 0.428e-6, 9.99e-12, 
	      6e-5, 3.5e-10, 30e-10, 0.4, 0.4);
    paramsets["MW3.1_60"] = p;

    p = -1.72, -0.322, 0.0254e-6, 0.438e-6, 3.20e-12, 
      5e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_50"] = p;

    p = -1.84, -0.132, 0.00898e-6, 0.489e-6, 2.90e-11, 
      4e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_40"] = p;

    p = -1.91, -0.125, 0.00837e-6, 0.499e-6, 4.15e-11, 
      3e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_30"] = p;

    p = -2.04, -0.111, 0.00828e-6, 0.543e-6, 5.57e-11, 
      2e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_20"] = p;

    p = -2.17, -0.0382, 0.00373e-6, 0.586e-6, 3.79e-10, 
      1e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_10"] = p;

    p = -2.25, -0.0648, 0.00745e-6, 0.606e-6, 9.94e-11, 
      0e-5, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["MW3.1_00"] = p;

    // for the LMC and SMC models, we adjust the abundance
    // normalization by 3.33 and 5.0, respectively, to account for the
    // fact that the WD01 models are normalized per H and we normalize
    // per mass of dust. Note that this correction is derived from
    // comparing the extinction values per H and per mass of dust in
    // the WD01 data files. It is inconsistent with the stated values
    // in the WD01 paper which says the abundance-limited values are
    // reduced by factors of 1.6 and 4.0 compared to the MW.
    p = -2.91, 0.895, 0.578e-6, 1.21e-6, 7.12e-17*3.33,
      0e-5*3.33, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["LMCavg_00"] = p;

    p = -2.99, 2.46, 0.098e-6, 0.641e-6, 3.51e-15*3.33,
      1e-5*3.33, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["LMCavg_10"] = p;

    p = 4.43, 0.0, 0.00322e-6, 0.285e-6, 9.57e-24*3.33,
      2e-5*3.33, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["LMCavg_20"] = p;

    p = -2.79, 1.12, 0.019e-6, 0.522e-6, 8.36e-14*5.0,
      0.0*4.0, 3.5e-10, 30e-10, 0.4, 0.4;
    paramsets["SMCbar"] = p;

    // The DL07 distributions have new PAH parameters and 0.92* the
    // number in the continuous distribution.
    p = -1.54, -0.165, 0.0107e-6, 0.428e-6, 9.99e-12*0.92, 
      6e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_60"] = p;

    p = -1.72, -0.322, 0.0254e-6, 0.438e-6, 3.20e-12*0.92, 
      5e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_50"] = p;

    p = -1.84, -0.132, 0.00898e-6, 0.489e-6, 2.90e-11*0.92, 
      4e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_40"] = p;

    p = -1.91, -0.125, 0.00837e-6, 0.499e-6, 4.15e-11*0.92, 
      3e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_30"] = p;

    p = -2.04, -0.111, 0.00828e-6, 0.543e-6, 5.57e-11*0.92, 
      2e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_20"] = p;

    p = -2.17, -0.0382, 0.00373e-6, 0.586e-6, 3.79e-10*0.92, 
      1e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_10"] = p;

    p = -2.25, -0.0648, 0.00745e-6, 0.606e-6, 9.94e-11*0.92, 
      0e-5, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_MW3.1_00"] = p;

    p = -2.91, 0.895, 0.578e-6, 1.21e-6, 7.12e-17*0.92*3.33,
      0e-5*3.33, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_LMCavg_00"] = p;

    p = -2.99, 2.46, 0.098e-6, 0.641e-6, 3.51e-15*0.92*3.33,
      1e-5*3.33, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_LMCavg_10"] = p;

    p = 4.43, 0.0, 0.00322e-6, 0.285e-6, 9.57e-24*0.92*3.33,
      2e-5*3.33, 4.0e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_LMCavg_20"] = p;

    p = -2.79, 1.12, 0.019e-6, 0.522e-6, 8.36e-14*0.92*5.0,
      0.0*4.0, 4e-10, 20e-10, 0.4, 0.55;
    paramsets["DL07_SMCbar"] = p;

  }
    
  // now it's easy, just set params
  std::map<std::string, params>::iterator i=paramsets.find(model);
  if (i==paramsets.end()) {
    std::cerr << "Error: Unknown WD01 model set: " << model << std::endl;
    throw false;
  }
  std::cout << "Using WD01 graphite size distribution " << model << std::endl;
  set_params(i->second);
};

void mcrx::wd01_silicate_distribution::set_params(const params& p) 
{
  alpha_ = p[0];
  beta_ = p[1];
  a_t_ = p[2];
  a_c_ = p[3];
  C_ = p[4];
};

mcrx::wd01_silicate_distribution::
wd01_silicate_distribution (const std::string& model,
			    const T_unit_map& units) : 
  wd01_distribution(units) 
{
  if (paramsets.empty()) {
    // initialize the paramsets
    // alpha, beta, a_t, a_c, C
    params p( -2.21, 0.300, 0.164e-6, 0.100e-6, 1.00e-13);
    paramsets["MW3.1_60"] = p;

    p = -2.10, -0.0407, 0.166e-6, 0.100e-6, 1.27e-13;
    paramsets["MW3.1_50"] = p;

    p = -2.10, -0.114, 0.169e-6, 0.100e-6, 1.26e-13;
    paramsets["MW3.1_40"] = p;

    p = -1.41, -11.5, 0.171e-6, 0.100e-6, 1.33e-12;
    paramsets["MW3.1_30"] = p;

    p = -1.43, -11.7, 0.173e-6, 0.100e-6, 1.27e-12;
    paramsets["MW3.1_20"] = p;

    p = -1.46, -10.3, 0.174e-6, 0.100e-6, 1.09e-12;
    paramsets["MW3.1_10"] = p;

    p = -1.48, -9.34, 0.172e-6, 0.100e-6, 1.02e-12;
    paramsets["MW3.1_00"] = p;

    // for the LMC and SMC models, we adjust the abundance
    // normalization by 3.33 and 5.0, respectively, to account for the
    // fact that the WD01 models are normalized per H and we normalize
    // per mass of dust.

    p = -2.45, 0.125, 0.191e-6, 0.100e-6, 1.84e-14*3.33;
    paramsets["LMCavg_00"] = p;

    p = -2.49, 0.345, 0.184e-6, 0.100e-6, 1.78e-14*3.33;
    paramsets["LMCavg_10"] = p;

    p = -2.70, 2.18, 0.198e-6, 0.100e-6, 7.29e-15*3.33;
    paramsets["LMCavg_20"] = p;

    p = -2.26, -3.46, 0.216e-6, 0.100e-6, 3.16e-14*5.0;
    paramsets["SMCbar"] = p;

    // the DL07 distributions are identical except 0.92*the number
    p = -2.21, 0.300, 0.164e-6, 0.100e-6, 1.00e-13*0.92;
    paramsets["DL07_MW3.1_60"] = p;

    p = -2.10, -0.0407, 0.166e-6, 0.100e-6, 1.27e-13*0.92;
    paramsets["DL07_MW3.1_50"] = p;

    p = -2.10, -0.114, 0.169e-6, 0.100e-6, 1.26e-13*0.92;
    paramsets["DL07_MW3.1_40"] = p;

    p = -1.41, -11.5, 0.171e-6, 0.100e-6, 1.33e-12*0.92;
    paramsets["DL07_MW3.1_30"] = p;

    p = -1.43, -11.7, 0.173e-6, 0.100e-6, 1.27e-12*0.92;
    paramsets["DL07_MW3.1_20"] = p;

    p = -1.46, -10.3, 0.174e-6, 0.100e-6, 1.09e-12*0.92;
    paramsets["DL07_MW3.1_10"] = p;

    p = -1.48, -9.34, 0.172e-6, 0.100e-6, 1.02e-12*0.92;
    paramsets["DL07_MW3.1_00"] = p;

    p = -2.45, 0.125, 0.191e-6, 0.100e-6, 1.84e-14*0.92*3.33;
    paramsets["DL07_LMCavg_00"] = p;

    p = -2.49, 0.345, 0.184e-6, 0.100e-6, 1.78e-14*0.92*3.33;
    paramsets["DL07_LMCavg_10"] = p;

    p = -2.70, 2.18, 0.198e-6, 0.100e-6, 7.29e-15*0.92*3.33;
    paramsets["DL07_LMCavg_20"] = p;

    p = -2.26, -3.46, 0.216e-6, 0.100e-6, 3.16e-14*0.92*5.0;
    paramsets["DL07_SMCbar"] = p;

  }
    
  // now it's easy, just set params
  std::map<std::string, params>::iterator i=paramsets.find(model);
  if (i==paramsets.end()) {
    std::cerr << "Error: Unknown WD01 model set: " << model << std::endl;
    throw false;
  }
  std::cout << "Using WD01 silicate size distribution " << model << std::endl;
  set_params(i->second);
};
